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High Frequency Asymptotics for 
the Spin- Weighted Spheroidal Equation 

Marc Casals* and Adrian C. Ottewill"^ 

Department of Mathematical Physics, University College Dublin, Belfield, Dublin 4, Ireland 

We fully determine a uniformly valid asymptotic behaviour for large auj and fixed 
m of the angular solutions and eigenvalues of the spin-weighted spheroidal differential 
equation. We fully complement the analytic work with a numerical study. 



I. INTRODUCTION 

By making use of the Newman- Penrose formalism, Teukolsky ( [1], [2]) showed that 
the equations describing linear scalar (spin-0), neutrino (spin-1/2) electromagnetic (spin-1), 
and gravitational (spin-2) perturbations of a general Type D background can be decoupled. 
Using the Kinnersley tetrad and Boyer-Lindquist co-ordinates, Teukolsky wrote the field 
equations in the Kerr background in compact form for the various spin fields, as one single 
'master' equation 
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where A = — 2Mr + and S = + cos^ 9. The field Qh and the source term Th 
are defined in [2]. The parameter h = 0, ±1/2, ±1,±2 refers to the helicity of the field. 
Following on Carter's work [3] for the scalar (spin-0) case, Teukolsky further showed that in 
the Kerr background the homogeneous decoupled equations can be solved by separation of 
variables: 
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The angular equation resulting from the separation of the Teukolsky equation is the so- 
called spin-weighted spheroidal differential equation and its regular solutions are the spin- 
weighted spheroidal harmonics (SWSH). In terms of a; = cos 6*, the spin- weighted spheroidal 
differential equation is 

d ( ,^ \ 2 2 im^hxf ^ ,2 



ax \ ax J 1 — x^ 

where hEimu denotes the eigenvalue and c = aw. It would be more logical to label the 
angular solutions and the eigenvalues by c rather than cu, but following convention we label 
them by uj. The eigenvalue for the case c = 0, corresponding to Schwarzschild space-time, 
is well-known to be 

hEi^^=^ = l{l + l)-h{h + l), (1.3) 

with regular solutions being the spin- weighted spherical harmonics [4]. 
The corresponding radial equation is 

A-^^ (a'^+^^^) - nVuRin.. = (1.4) 

where the potential is given by 

2ih(r-M)K-K^ , , , , 

hV = ^ ^ Aihur + h\mc. (1.5) 

with K = {r'^ + a^)a; — am. The separation constants are related by 

h>^imoj = hEimoj - h{h 1) + - 2mc (1.6) 

The differential equation (1.2) has two regular singular points at x = ±1 and one essential 
singularity at x = oo. We are only interested in solutions for real values of the independent 
variable 9 corresponding to the interval x G [—1, +1]. We henceforth restrict x to this range 
and therefore we have only to consider the two regular singular points at x = ±1. The 
differential equation (1.2), together with the boundary condition that its solution hSimuj{x) 
is regular for x G [—1,-1-1], defines a parametric eigenvalue problem, with parameters c, m 
and h. The physical requirements of single- valuedness and of regularity at x = ±1 requires 
that / and m are integers with |m| < /. Stewart [5] showed that the SWSH form a strongly 
complete set if c is real while he could only prove weak completeness if c is complex. 

In this paper we study the asymptotic behaviour for high frequency of the solution and 
eigenvalues of the spin-weighted spheroidal differential equation. Following standard con- 
ventions, we refer to 'high frequency' in relation to the angular solution and eigenvalues 
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when in fact what it is meant is large c(= acj), where a is the angular momentum per unit 
mass of the rotating black hole and u is the frequency of the mode. 

The high frequency approximation of the spin-weighted spheroidal equation is a partic- 
ularly important subject that has been left unresolved thus far, except for the spin-0 case, 
due to its difficulty. This asymptotic study is important when considering both classical and 
quantum perturbations. In the classical case it is important, for example, when calculating 
gravitational radiation emitted by a particle near the black hole since the typical time-scale 
of the motion is short compared to the scale set by the curvature of the black hole. In the 
quantum case its importance lies in the fact that the high frequency limit is at the root of 
the divergences that the expectation value of the stress-energy tensor possesses. The correct 
subtraction of the divergent terms from the expectation value of the stress-energy tensor 
is extremely troublesome in curved space-time, particularly in one that is not spherically 
symmetric. As the divergent terms arise from the high frequency behaviour of the field, 
knowledge of this behaviour is fundamental in such a subtraction. This limit has also been 
recently considered in the Kerr background in the context of quasinormal modes (see [6]). 

All analysis in this paper has been performed for general spin, so that it applies to the 
scalar, neutrino, electromagnetic and, in particular, gravitational perturbations, which are 
of great interest in astrophysics. However, we should note that the asymptotic study in this 
paper is valid for fixed m as c tends to infinity, a fuller understanding of the asymptotic 
behaviour of the solution would require an anlysis uniform in m. 

In the remainder of this introductory section we discuss the results for high frequency 
asymptotics of SWSH that have been obtained in the literature up until now, show their 
shortcomings and outline what our new results achieve. In the next section we lay down the 
basic theory that we use in the following sections. In Sections III A, III B, III C and IV we 
fully determine the aymptotic behaviour of the angular solution that is uniform in x and the 
asymptotic behaviour of the eigenvalue. In Section V we describe the numerical method and 
programs used to obtain the numerical results, which in the last section we show, analyze 
and compare to our asymptotic results and to numerical results in the literature. 

Different authors have obtained high-frequency approximations to the solution and eigen- 
values of the spheroidal differential equation, which results from the spin- weighted spheroidal 
differential equation when /i = 0. Erdelyi et al. [7], Flammer [8] and Meixner and Schafke 
[9] have all done so using the fact that the spheroidal differential equation becomes the 
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Laguerre differential equation in the high-frequency hmit. 

Breuer [10] was the first author to study the high-frequency behaviour of the spin- 
weighted spheroidal harmonics. Based on the work on the spin-0 case by the above authors 
he related the solution of a transformation of the spin-weighted spheroidal equation for large 
c and finite m to generalized Laguerre polynomials. His work, however, was fundamentally 
flawed as it assumed that the solution was either symmetric or antisymmetric under 9 — >■ 
TT — 6, which is only true for spin-0. 

Breuer, Ryan and Waller [11] (hereafter referred to as BRW) corrected this error and 
further developed this study by first relating the SWSH to the confluent hypergeometric 
functions and then reducing them to the generalized Laguerre polynomials by imposing 
regularity far from the boundary points x = ±1, where x = cos 6. Unfortunately, their 
study of the high-frequency behaviour was also flawed and incomplete. The behaviour 
for high frequency of both the spherical functions and the eigenvalues obtained by BRW 
depend critically on a certain parameter 7 (called a in that paper) which they were unable 
to determine for the case of non-zero spin. 

BRW obtained the analytic value of 7 for the spin-0 case, however for non-zero spin 
they could only calculate it numerically for a handful of sets of values of {/, m} for spin- 
2. BRW achieved this numerical calculation for the spin-2 case by matching the high- 
frequency asymptotic expression for the eigenvalue that they obtained with the expression 
for the eigenvalue given by Press and Teukolsky [12] valid for low frequency. Not only their 
analytic expressions for both the spherical solution and the eigenvalue for high frequency 
were thus left undetermined, but also their expressions for the spherical solution are only 
valid sufficiently close to the boundary points x = 1 and x = —1, not for the region in- 
between them. This results in the possibility that a zero of the solution near x = 0, away 
from X = ±1, be overlooked. Furthermore, and crucially, their assumption that the confluent 
hypergeometric functions should reduce to the generalized Laguerre polynomials by imposing 
regularity far from the boundary points is not correct. The reason why it is not correct is 
that in the cases for which the confluent hypergeometric function diverges far from one of the 
boundaries, the coefficient in front of it decreases exponentially with c so that the solution 
remains finite in the whole region x G [—1, +1]. We believe that the reason why they were 
not able to analytically determine the value of the parameter 7 is because they ignored the 
behaviour of the solution far from the boundaries, thus overlooking a possible zero, and 
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wrongly imposed regularity. 

The study of the behaviour of the solution and eigenvalues of the spin- weighted spheroidal 
equation for high frequency and finite m has not been developed any further by these or any 
other authors and therefore BRW's work is where this study stood until the present paper. 

In this paper we correct and complete BRW's study for high frequency and finite m. 
We thus obtain an asymptotic solution for large frequency to the spin-weighted spheroidal 
equation which is uniformly valid everywhere within the range x G [—1, +1], not just near 
the boundaries. We also analyze the existence and location of a possible zero of the solution 
near x = 0. We analytically determine the value of 7 by matching the number of zeros that 
our asymptotic solution has with the number of zeros the SWSH has. As a consequence, the 
asymptotics of the eigenvalue in the same limit also become fully determined. Finally, we 
have complemented all the analytic work with graphs produced with data that we obtained 
numerically. The graphs show the behaviour of the eigenvalues for large frequency and how 
they match with Press and Teukolsky's approximation for low frequency. They also show 
the behaviour of the SWSH in this limit and the location of its zeros. 

II. SYMMETRIES OF THE SPIN- WEIGHTED SPHEROIDAL DIFFERENTIAL 

EQUATION 

Certain symmetries of the spin-weighted spheroidal equation (1.2) are immediate: the 
equation remains invariant under the change in sign of two quantities among [s, {m,uj),x], 
where we are considering that (m, cu) constitutes one single quantity, i.e., m and u change 
sign simultaneously. As a consequence, the SWSH satisfy the following symmetries, where 
the choice of signs ensures consistency with the Teukolsky-Starobinskii identities below: 

hSinUO) = (-l)'+'"-,5z™^(7r - 9) (2.1a) 

hSlmUO) = {-ly'-^.Sl.rn-U^ " 0) (2.1b) 

hSimM = (-l)''+"'_,.Si_„_^(e) (2.1c) 
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Here any one symmetry follows from the other two. The eigenvalues must consequently also 
satisfy the symmetries: 

hElmuj = -hElmoj (2.2a) 
hElmuj = hEl^m.-ui (2.2b) 

The SWSH with helicity h is related to the SWSH with helicity —h via the Teukolsky- 
Starobinskii identities [13]. We start by defining the operator 

/ \ rn 

CV^=deTcsm9± - + ncote. (2.3) 

smO 

Then for spin-i, the Teukolsky-Starobinskii identities may be written as 

^l+lSimuj = ~lBlmu}^lSlmuj (2.4a) 
2^2 2 2 



— ^^ImiU + ^Simuj (2.4b) 

where 



^^Imui — -l^lniLV (2.5) 

For the spin-1 case, the Teukolsky-Starobinskii identities may be written as 

^O^l+lSlmuj = llSlmuj-lSlmi^ (2.6a) 

Sim. (2.6b) 

where 

A. = -lAL. + 4mc - 4c' (2.7) 
Finally for the spin-2 case, the Teukolsky-Starobinskii identities may be written as 

C-lCoCiC2+2Slmuj = 2^lmu)-2Slmu) (2.8a) 

C^_lC\C\C\-.2Slmu) = 2^lmuj+2Slmu} (2.8b) 



where 



^"imui — -2^'lmuji-2^lmui + 2)^ — SC^ -2hmui | (^1 [5-2^lmuj + 6] — 12| 



(2.9) 



The signs of ^Bimuj, iBimu and 2'Szmw are arbitrary, but we will take them to be both positive. 
With this convention, (2.4), (2.6) and (2.8) agree with the sign in the symmetry (2.1a) of 
the angular function. 
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III. THE UNIFORM ASYMPTOTIC SOLUTION 

In the rest of this paper we follow the approach to boundary layer theory as presented 
by Bender and Orszag [14]. The asymptotic solution that is a valid approximation to the 
solution of the differential equation from the boundary point ±1 until x ~ ±l + 0(c'^), where 
— 1 < (5 < 0, is called the inner solution. The region within which an inner solution is valid 
is a boundary layer. As we shall see, for the large frequency approximation of the spin- 
weighted spheroidal equation, there are two boundary layers within the region x G [—1, 1], 
one close to x = —1 and one close to x = +1. Close to the boundary points the SWSH 
oscillate rather quickly in x, and indeed it is there where all the zeros of the function are 
located (with the possible exception of one). 

The asymptotic solution that is a valid approximation to the solution of the differential 
equation in the range — 1 + 0(c~^) <^ x ^ +1 — 0(c~^), is called the outer solution. This 
range comprises not only the region in between the two boundary layers but also a certain 
region of both boundary layers. This region where both an inner solution and the outer 
solution are valid is the overlap region, and it is there that the outer and inner solutions are 
matched. 

We shall see that in between the two boundary layers the function behaves rather 
smoothly, like a coshx or a sinhx, so that the SWSH may have at the most one zero, 
which will turn out to lie close to x = 0. The behaviour of the outer solution is important 
despite its smoothness because when matching it with the inner solutions it will allow us to 
find an asymptotic solution which is uniformly valid throughout the whole range of x. The 
outer solution is also necessary in order to find out whether or not the uniform solution has 
a zero close to x = and, if it does, to calculate the analytic location of the zero. 

This is a key feature that singles out the scalar case from the others: for the spin-0 case 
the differential equation (1.2) is clearly symmetric under {x ^ — x} and therefore, depending 
on its parity, it will have a zero at x = or not. On the other hand, for the case of spin non- 
zero, the differential equation does not satisfy this symmetry but it does remain unchanged 
under the transformation {x ^ — x, h ^ —h} instead. There is therefore no apparent reason 
why it should have a zero near the origin. The outer solution is important for the case of 
spin non-zero and not for spin-0 since, as we shall see, the differential equation that the 
outer solution satisfies is symmetric under {x <-> — x, h ^ —h} to leading order in c. 
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As noted above, the differential equation (1.2) has singular points at x = ±1. By using 
the Frobenius method it can be found that the solution that is regular at both boundary 
points a; = +1 and —1 is given by 

hSim.{x) = (1 - + xfhyi^^{x) (3.1) 

where 

\m + h\ ^ \m-h\ „x 

a = , p = (j.^j 

2 2 

and the function hVimujix) satisfies the differential equation 

l-x'')^-2[a-p + {a + P + l)x] ^+ 
ax'^ ax 

+ hEimoj - {a + (3){a + (3 + l) + (?x^ - 2/icx|/j?/im^(a;) = 

(3.3) 

A. Inner solutions 

BRW obtained an expression for the inner solution for general spin in terms of an unde- 
termined parameter 7. In this section we summarize and present their results in a compact 
way. 

By making the variable substitution u = 2c(l — x), equation (3.3) becomes 
u — — h (2tt + 1) 



du"^ du 

- ^ |m + 2/l - ^ [c^ - (a + /3)(a + /3 + 1) + hElrrJ\ | hVlnu^- (3.4) 

- ^ |n + 2(a + /5 + 1)^ - ^-n + hu^ j = 

It is clear from this equation that the leading order behaviour of hEimtj for large c must be: 

hEimu, = -c^ + -fc + 0{l). (3.5) 

If its leading order were not — c^, there would then be a leading order term +^Chyimuj in the 
equation that it could not be matched with any other term. Lower order terms for hEimuj 
are given in BRW. It is crucial to know the value of the parameter 7, as it determines how 
the angular function behaves asymptotically to leading order in c. At this stage, 7 is an 
undetermined real number; we will determine its value later on. 
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Using the asymptotic behaviour (3.5) and letting c — > oo, the terms in (3.4) of order 
0(c~^) can be ignored with respect to the other ones and, to leading order in c, the function 
hVimw satisfies 

^~dv? ^ '~(hl 4 + - 7) hVlrnu: = 0. (3.6) 

The solution of this differential equation that satisfies the boundary condition of regularity 
at X = +1 is related to the confluent hypergeometric function: 

hy'rZ''' = hCim.e-'''\Fi ((\m + h\ + h + l)/2 - 7/4, \m + h\ + 1, n) (3.7) 

where hCimu is a constant of integration. 

Similarly, if we instead make a change of variable u* = 2c(l + x) in equation (3.3), due 
to the {x ^ —X, h ^ —h} symmetry we obtain 

/^^/;^:"' = hD,m^e-^*/\F^ ((|m - h\ - h + l)/2 - 7/4, \m - h\ + 1, n*) (3.8) 

as the solution that is regular at x = — 1. 

We use the following obvious notation to refer to the solutions of the spin-weighted 
spheroidal equation that correspond to the inner solutions of (3.6): 

The inner solution hSl^^^ is only a valid approximation in the region from the boundary 
point ±1 until ±1 — x ~ 0{c^) with — 1 < 5 < 0. The reason is that in the step from (3.4) to 
(3.6) we have ignored terms with u^*^ /c with respect to terms of order 0(1), and therefore 
the inner solution has been found for ±1 — x~n'f*J'/c^ 0(1) and so we must have 6 < 0. 
On the other hand, we are not ignoring u with respect to the 0(1) term (2/1 — 7) iii equation 
(3.6), so that it must be m ~ 0{c^'^^) with 6+1 >0. From the fact that we are not ignoring 
{2h — 7) with respect to u it does not follow that 6 + 1 < 0, since the inner solution is valid 
at the boundary point x = +1, where u = 0. That is, the term {2h — 7) cannot be ignored 
with respect to u for all x from +1 up to +1 — x ~ 0{c^), even if 5 + 1 > 0. A similar 
reasoning applies to u*. 

We therefore have one boundary layer comprising the region in x from —1 to (— 1 — x) ~ 
0{c^) and another boundary layer from (+1 — x) ~ 0{c^) to +1. 
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To leading order in c the solution to the spin-weighted spheroidal equation which is valid 
within the two boundary layers is given by 

.«,-^(i-.r(i..)''r'^'""""''"''""""'"' (3.9) 

/xAm.e-"*/\Fi(-p',2/3 + 1,M*) X < 



where we have defined 



-(|m + /i| + /i + l)/2 + 7/4 
-(|m-/i| -/i + l)/2 + 7/4 



(3.10) 



BRW then require that p^p' E Z+ in order that the inner solution hSl^^ is regular at 
X = 0, where u,u* — * oo. Correspondingly, they replace the confiuent hypergeometric 
functions iFi{a,b,x) by the generalized Laguerre polynomials L^^~^\x). As we shall see, 
this is erroneous: p, p' G Tj^ is not a necessary condition for regularity since in the cases for 
which this condition is not satisfied, the coefficients hCimw and hDimuj diminish exponentially 
for large c in such a way that hS}'^^ remains regular. 



B. Outer solution 

We now proceed to find the outer solution of the spin-weighted spheroidal differential 
equation. We first make the variable substitution 

y{x)=g{x)exp / ^^^±i^±l±ll^rfx = (7(x)(l - x)-(2-+l)/2(l + 3;)-{2/3+l)/2 (g^^^ 

J 1 — x^ 

which transforms equation (1.2) into 

g"{x) + f{x,c)g{x) = Q (3.12) 

where 

G(x,c) (a + (3 + l)(l-x'^) + 2x[a-(3 + (a + (3 + l)x] 
j[x,c) = + 



1 — (1 — x^)^ 

[a- P + {a + P + l)x] 



(3.13) 



(1 -X2)2 

and G{x,c) is the coefficient of hUimuj in (3.3), i.e. 



G{x, c) = hEimoj - {a + p){a + p + 1) + c^x^ - 2hcx 



(3.14) 
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We now perform a WKB-type expansion: g{x) = e^^^\ This change of variable converts 
equation (3.12) into 

g"{x) + g'ixf + f{x,c) = (3.15) 
Performing an asymptotic expansion of f{x, c) in c we find 

fix, c) = foixy + h{x)c + 0(1), (3.16) 

with 



/o(x) = -1, h{x) = "7 ;;^ (3.17) 



2{q - hx) 
1 -x2 

where we have used the asymptotic expansion of hEimui in c and we have also introduced the 
parameter q = 7/2. We will prove in Section IV that q must be an integer. It is clear that 
to leading order in c the outer solution is symmetric under {x ^ — x}. We are avoiding any 
possible turning points by assuming that f{x, c) 7^ for x values of interest. This condition 
is clearly satisfied if c is large enough. 

Next we perform an asymptotic expansion of Q{x) in c. We do not know a priori what 
the leading order is, and so we will determine it with the method of dominant balance. Let 
the expansion of Q{x) for large c be Q{x) = ho{c)Qo{x) + o(/io(c)). On substituting the 
asymptotic expansions for f{x,c) and Q{x) into (3.15) we obtain 

ho{c)g'^{x) + hoicf {g'o{x)f + cVo(x) + o{ho{cf) + o{c^) = (3.18) 

We could try and cancel out the c^/o(x) term with /;,o(c)^o (^)i ^^at would give ho = c^, but 
then /io(c)^g (x) would be subdominant to /io(^o)^- "^^^ other option is to cancel the c^fo{x) 
term with /io(^o)^ instead. This gives Hq = c, which works. We therefore have that 

G{x) = cQoix) + g^ix) + O(c-i) (3.19) 

The resulting equation for the leading order term in g is 

[g'^(x)]' + Ux) = 0, (3.20) 

the solution of which is go = ±{x — Xq) . The equation for the next order in c is 

g'^ix) + 2g',{x)g[{x) + h{x) = o, (3.21) 

which gives 



Qi = T 



log(l + x) — log(l - x) 



(3.22) 
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The physical optics approximation for the outer solution is therefore given by 

HyZ.{^) = (1 - xr^'-^'^l\l + a:)-(^^+^)/^x 

(3.23) 

where the constant Xq has been absorbed within h^imuj and hBimu- This solution is valid in 
the region -1 + 0(c"^) < x < +1 - 0(c"^). 

C. Matching the solutions 

We have found three different solutions. One of the two inner solutions is valid in the 
region — 1 < a; < — 1 + 0{c^) for any 5 such that — 1 < 5 < 0, and the other one for 
+ 1 - 0{c^) < X < +1. The outer solution is valid for -1 + 0(c~^) < a; < +1 - 0(c~^). 
Clearly all three solutions together span the whole region — 1 < x < +1. There are also two 
regions of overlap, one close to -1 and one close to +1, where both the outer solution and 
one of the inner solutions are valid. We can proceed to match the solutions in these regions 
and we will do so only to leading order in c as matching to lower orders would not bring any 
more insight into the behaviour of the SWSH. When the matching is completed to leading 
order, the two overlap regions are given one by 0(c~^) ^ 1 + x < 0{c^) and the other one 
by 0(c"^) ^ 1 — X < 0{c^). For the overlap regions to exist it is therefore required that we 
choose a 5 satisfying —1 < 5 < 0. 

In order to obtain an expression for the inner solution in the overlap region, we expand the 
inner solution for u^u* ~ oo. For that, we need to know how the confluent hypergeometric 
functions behave when the independent variable is large. From [15] we have 

when z = \z\e^^ with — 7r/2 < < Sir /2, which includes the case we are considering: = 0. 
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This means that the inner solution vahd close to x = +1 behaves like 

r(|m + h\ + l) [2c(l - e+^(^-") 

inn.+l f-y J ^{~P) 



p^Z+U {0} 



inn, "hi ^ 



r{\m + h\ + l)e— P [2c(l - x)r e^^^ p e Z+ U {0} 



T{\m + h\ + l+p) 

, +oo) 

(3.25) 

The behaviour of the inner solution valid close to x = — 1 is similarly obtained by simulta- 
neously replacing x with — x, h with —h (which also implies replacing p by p') and hCimuj 
with hDimuj above. 

On the other hand, in order to obtain an expression for hV'imuj valid in the overlap region 
we perform a Taylor series expansion around x = +1 or — 1 depending on where we are 
doing the matching, and keep only the first order in the series: 

a) Around x = +1. 

To first order in (1 — x): 

kVZ'M - .A_(l -x)[+(<^"^"i)/2""]2[-('^+^+i)/2-^]e+-+ 

(3.26) 

+ f^BimM - x)[-('?"^+l)/2-"l2[+(«+'^-l)/2-;3]g-c. ^ 

By matching the inner and outer solution in the overlap region 0(c~^) <C 1 — x < 
O(c^), i.e., by matching equations (3.25) and (3.26), we obtain the following relations 
depending on the value of p: 

al) if p ^ Z+ U {0}: 

,B,^ = 2hto+^-i)/^+/^l ^('";+ ^1 + ^^ (2c)^^-l».+^l-ile+c,q_ ^^-^^^ 
a2) iipeZ+U {0}: 

= 2[+(^+'^+^)/2+^] rn ^'"^ti^itl^ , e-^(2c)^e-,q^. (3.28) 
r(|m + h\ + 1 + p) 

b) Around x = — 1 (similar to the x = +1 case). 
To first order in (1 + x): 

ky?l.{x) - .A,_(l + x)["(^+'^+^)/^-'^]2[+(^-'^-i)/^""]e+^^+ 



(3.29) 
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S 



inn, — 1 



^inn,+l 



S' 



out 



:'match, — 1 



^match,+l 



— 1 \> fc> 'o^ 

OVJ 0^ 0^ 







fc-^ \> +1 
^ x^ 



FIG. 1: Regions of validity in the x axis of the various approximations to the SWSH for large c. It 
must be— l<e<5<0. For clarity, the mode labels have been dropped, ^^^tchjii pgfgpg f^j^g 
asymptotic approximation valid in the overlap region (red) close to x = ±1. The uniform solution 
is constructed as 5"°'^ = 5°"* + 5"^°'+^ + S''^'^~^ - S-^atch.+i _ ^-match-i^ 



bl) if ^ Z+ U {0}: 

hBlmuj = 



(3.30) 



-p' 

b2) if G Z+ U {0}: 

= 2[+('^-^+^)/^+-l - ^ + e-^'(2c)^'e-.An.. (3.31) 

1 (|m — Ai| + 1 + p') 

From the above matching equations we can obtain a uniform asymptotic approximation 
to hSirrnu valid throughout the whole region x G [—1, +1] and also find out where the zeros 
of the function are. The uniform asymptotic approximation is obtained by adding the outer 
and the two inner solutions, and then subtracting the asymptotic approximations in the two 
overlap regions since these have been included twice. Figure 1 depicts the region of validity 
of the various asymptotic solutions for large c that we have obtained. 

We can distinguish three cases: 

p,p' ^Z+U {0} 

From equations (3.27) and (3.30) it must be hAmuj = hBimw = = hCimw = hDinnj, so 
this case is the trivial solution and we discard it. 
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p e Z+ U {0} and p' ^ Z+ U {0} , or vice-versa 

Either hA^imu or hBimui is equal to zero (but not both), so that the function hSimu cannot 
have a zero close to x = 0. All the zeros, if there are any, of hSimw are zeros of the inner 
solutions and thus they are located inside the boundary layers, close to x = ±1. 



In this case we can already directly obtain the uniform asymptotic approximation, 
up to an overall normalization constant hCimuj'- 



x) + 



hSZl = hQmM - xni + x)^|e-^(i-^)iFi ( - p, 2a + 1, 2c(l 

r{2a + 1) T{-p') .^^p, ,p^p/^2/;+l -2co(q+/3-a) -c(l+x) 

+ r(2a + l+p)r(2/5 + l)' ^^'^ ' ^ ' 

X iFi f - p', 2/5 + 1, 2c(l + x)) + 2[+(''+'*+i)/2+/31 r(2a + l) ^ 
V / r(2a + 1 + 

_2[+to-/^-i)/2-a]^l + x)-(''+'^+^)/2-/3j I ^hen p e Z+ U {0} and ^ Z+ U {0} 

(3.32) 

The uniform approximation when p ^ Z"*" U {0} and p' G Z'^ U {0} may be obtained by 
making the substitutions x ^ —x and h ^ —h (which imply the substitutions a ^ l3 and 
p «-> p') in (3.32). 

The irregularity arising from e~^^^~^^\Fi{—p' ,2P + l,2c(l + x)) ~ e^'^ (ignoring factors 
independent of x and c) in the limit x ^ +1 and c ^ +oo prompted BRW to discard the 
case p' ^ Tj^ U {0}. It is clear from (3.32), however, that this irregularity is nullified by the 
factor e~'^'^ in front of it, brought in by the coefficient hDimu)- Note that despite the factor 
e"^"^, close to X = — 1 this term (which is part of the inner solution valid in the boundary layer 
there) is not dominated by the first term in (3.32) (which is the inner solution valid in the 
boundary layer near x = +1). The reason is that e~'^^^~^^\Fi{—p' , 2(3 + 1, 2c(l + x)) ~ e~'^^ 
and e~^'^^~^\Fi{—p^ 2a + l, 2c(l— x)) ~ e""^^ where both limits are x ^ — 1 and c — > +oo and 
we have ignored factors independent of x and c. In the boundary layer around x = ±1, the 
asymptotic approximation valid in the overlap region close to x = =f1 cancels out the inner 
solution hS^^^^^ in expression (3.32). Similarly, in the same boundary layer, the asymptotic 
approximation valid in the overlap region close to x = ±1 cancels out the outer solution. 
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so that only hSl^^^^ contributes to the uniform approximation in that boundary layer. A 
similar reasoning can be applied to the case p ^ U {0}. 

p,p'ez+u{o} 



In this case, apart from the overall normalization constant there is another unknown 
constant. We are going to determine this extra unknown by imposing the appropriate 
parity under {x ^ —x,h ^ —h}. Using the Teukolsky-Starobinskii identities (2.4), (2.6) 
and (2.8) together with the symmetry (2.1a) in the outer solution (3.23) we obtain 



2 



2 _ /_-|^\(/+m) 72 



D 



Imuj 



\{l+m) 



V2^ 



m 



for spin- 1/2, 

— iDlmijj 



-iQ 



Imuj 



m + 1/2 
m-1/2 1 

\^^q - m a/c 



when m> +\ 



when m < 



(3.33) 



'2^{q~m- l)(g-m + 1) 



(-1) 



{l+m) 



for spin-1 and 



. 2y/(g -m-l){q-m + l)c 



Y-^{l+m) +2Clmuj 
-2Clrau> 



when m > +1 
when m = 
when m < —1 



(3.34) 



4\/ {q — m — l)(g — m + l)(g — m — 3)(g — m + 3) 



(m + 2)(m + l)m(m 
Vg(g-2)(g-4) 



3v/^T2 
V{q-3)(q-l) 

V(g + 3)(g + l) 
3v/F^ 



v/g(g + 2)(g + 4)c 

(m + l)m(m — l)(m — 2) 

4\/ {q — m — l){q — m + l){q — m — 3){q — m + 3) 



when m > +2 

when m = +1 

when m = 

when m = — 1 

when m > —2 

(3.35) 
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for spin-2. Equations (3.33)-(3.35) have been obtained without imposing any restrictions on 
the values of p or p' and might therefore seem to contradict the resuh from (3.28) and (3.30) 
[or (3.27) and (3.31)] giving an exponential behaviour with c for the ratio hDimuj/ hCimuj for 
the case p G Z+ U {0} and p' ^ Z+ U {0} [or viceversa]. We shall see in the next section, 
however, that equations (3.33)-(3.35) can only actually be applied to the casep,p' G Z+U{0} 
so that there is no such contradiction. 

We can already determine in what cases the outer solution has a zero. Clearly, from 
equations (3.28), (3.31) and (3.33)-(3.35), the ratio between the coefficients h^imuj and hBimuj 
is proportional to a power of c, where the constant of proportionality does not depend on 
c. It then follows from the form (3.23) of the outer solution that one exponential term will 
dominate for positive x and the other exponential term will dominate for negative x, when 
c oo. Therefore the outer solution does not possess a zero far from a; = for large c. The 
outer solution has a zero if h^imui and hBimu have different sign and it does not have a zero 
otherwise. From equations (3.10), (3.28), (3.31) and (3.33)-(3.35) we have: 



Furthermore, we can calculate what the location of the zero of the outer solution is to 
leading order in c: by setting the outer solution (3.23) equal to zero and using (3.28) and 
(3.31) (since we have already seen that if p and/or p' ^ Z"^ U {0} the outer solution does not 
have a zero) we obtain that for large frequency the zero is located at the following value of 
x: 



Clearly, there is one zero in the region between the two boundary layers tending to the 
location x = as c becomes large if hAimu and hBimu have different sign and there is not a 
zero if they have the same sign. 




(3.36) 




i~h+a~f3) r(|m -h\ + l)T{\m + h\ + 1 + p) ^^^(p/_p) , . ^p'_p) hPlmu^ 

r(|m + ;i| + l)r(|m-/i| + l+pO ^ ^ hCi^^ 



) 

3.37) 
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Finally, the uniform asymptotic approximation for this case is: 



1 + a:)^|e-'=(^-^\Fi - 2a + 1, 2c(l - + 
+ ^^e-^(^+^)iFi ( - p', 2/5 + 1, 2c(l + x)) + 

^l(q+h+l)/2+(3] r(2a + 1) ^^^p. ,p ^ 

r(2a + l+p) ^ ^ 

X [(1 - 3;) + (g-?»-l)/2-a(^]^ _^ ^^|-(g+h+l)/2-/3 _ 2-[{g+h+l)/2+/3] (^-|_ _ ^y{q-h-l)/2-a^ _^ (3.38) 
/xA™^^[(g_/i+l)/2+Q] r(2/3 + 1) i^p,,^.p> 

^ nCimJ r(2/5 + l+p0 ^'"^ 

when p,p' G Z+ U {0} 

where the ratio between hDimw and hCimw is given by (3.33)-(3.35). 

Similar cancellations to the ones for the case p G Z+ U {0} and p' ^ Z+ U {0} occur in 
the present case for the uniform solution (3.38). The only difference is that now, in the 
boundary layer around x = ±1, the asymptotic approximation valid in the overlap region 
around x = ^1 only cancels out part of the outer solution. The other part of the outer 
solution, however, is exponentially negligible with respect to the inner solution hSl'^^^. 



IV. CALCULATION OF 7 

In order to finally determine the value of 7 we only need to impose that our asymptotic 
solution must have the correct number of zeros. BRW give the number of zeros of the SWSH 
for non-negative m and h. Straightforwardly generalizing their result for all possible values 
of m and h using the symmetries of the differential equation, we have that the number of 
zeros of hSimuj is independent of c and for x G (—1, 1) is equal to 



(4.1) 



The number of zeros of the confluent hypergeometric function is also needed, and that is 
given by Buchholz [16]: 





m 


for m 




\h\ 


for \m\ 
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The number of positive, real zeros of a, b, z) when 6 > is 

— [— al for + oo > a > , ^ 

\ - (4.2) 

for > a > — oo 

where [n] means the largest integer < n. 

Since the confluent hypergeometric functions are part of the inner solutions and the 
region of validity of these solutions becomes tighter to the boundary points as c increases, 
the zeros of iFi{—p,2a + l,n) are grouped together close to x = +1, and likewise for 
iFi(— p',2/5 + 1,M*) close to x = —1. Apart from these zeros, for large c the function hSimu 
may only have other zeros ai x = ±1 and/or at a; = Xq. The possible one at x = Xq is 
not due to the confluent hypergeometric functions but to the outer solution. We define the 
variable zq so that it has value +1 if hSimuj has a zero at x = xq and value if it does not. 

From equation (3.10) we see that p' = p + (|m + h\ + 2h — \m — h\)/2, and therefore if 
either p or p' is integer then the other one must be integer as well. But, as we saw in Section 
III C, at least one of p and p' (if not both) must be a positive integer or zero. Therefore 
both p and p' must be integers and at least one of them is positive or zero. It also follows 
from (3.10) that 

^ = 2{p + p') + 2+\m + h\ + \m-h\=2q (4.3) 

where it is now clear that g G Z. 

Requiring that the number of zeros of the asymptotic solution coincides with the number 
of zeros of the SWSH results in the condition 

-(|m + /i| + /i + l)/2 + g/2 for q>\m + h\ + h+l 
foT q < \m + h\ + h + 1 

^ -{\m- h\- h + l)/2 + q/2 ioT q>\m-h\-h + l 
for g < |m — /i| — /i + 1 



I — \m\ for |m| > \h\ 
I — \h\ for |m| < 

From (4.4) and the fact that zq = when either p or p' ^ Z+ U {0} as seen in Section 
III C, we obtain the value of q in all different cases: 
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/ - |m| for \m\>\h\ I (^\rn + h\ + \m - h\) 

Q = \ > H 7: - zo 

I - \h\ for \m\ <\h\ j 2 ^45g^) 

if />/i,/2 (i.e.,p,p'GZ+U{0}) 

{/ — \m\ for \m\ > \h\ 
'~''} + \m + h\ + h + l 
I - \h\ for |m| < \h\ J (4.5b) 

if l<l2 (i.e.,pG,p'^Z+U{0}) 

{/ — \m\ for Iml > 
> + |m - /i| - /i + 1 
/ - \h\ for |m| < \h\ J (4.5c) 

if l<h (i.e., p G Z+ U {0}) 

where 

li={ !""!",''",' } +i\m + h\-\m-h\)/2 + h 

+ {\m-h\ - \m + h\)/2-h 

By requiring in (4.5a) that q must also satisfy (3.10) and bearing in mind that zq can 
only have the values or 1, it must be 

{0 for / — li even 
(4.6) 
1 for l-k odd 

where I2 instead of li could have been used, since one is equal to the other one plus an even 
number. 

It can be trivially seen that if li has an allowed value, i.e.. 



m 


for 


m 


>\h\ 


\h\ 


for 


m 


<\h\ 


m 


for 


m 


>\h\ 


\h\ 


for 


m 


<\h\ 



1 


for 


|m| > \h\ [ 


( \h\ 


for 


m < \h\ J 



h>{ >, (4.7) 



then I2 does not, and vice-versa, so that cases (4.5b) and (4.5c) are mutually exclusive. 

Clearly, when I < li or I < I2, for fixed h and m, as / is increased by 1 the corresponding 
value of q is also increased by 1, so that two different values of / correspond to two different 
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values of q. However, once the threshold / > max(/i,/2) is reached, every increase of 2 in 
/ will involve the subtraction of an extra 1 in (4.5a) via Zq, so that its corresponding value 
of q will be the same as for the previous /. Therefore, in the region / > max(/i,/2), every 
value of q will correspond to two consecutive, different /'s: the two corresponding SWSH's 
will have the same number of zeros and behaviour close to the boundary points, but one 
will have a zero at x = Xq and the other one will not. 

Another feature that can be seen is that, for h = ±i, the case / < Zi or I2 (i.e., p or 
p' ^ Z+ U {0}) implies g — m = Oor/<|m| when m > | and m < — | respectively, so that 
(3.33) is not apphcable to these already mentioned in the previous section. 

Similarly, for h = ±1, the case / < /i or I2 implies q — m = ±1 org = +1 when m > 1 
and m = respectively, so that (3.34) is not valid for these cases. When m < — 1 it follows 
from (4.5b) and (4.5c) that / < /i or I2 requires / < \m\, which is not allowed. 

For h = ±2, / < /i or I2 implies 



m - 


-q = ±l,±3 


when 


m 


> 2 


Q = 


0,2,4 


when 


m 


= 1 


Q = 


1,3 


when 


m 


= 


Q = 


2 


when 


m 


= -1 


I < 


m 


when 


m 


< -2 



V. 



so that (3.35) is not valid then. 

Note that the scalar case is obtained from our formulae as a particular case. Setting 
/i = in the equations above we have li = I2 = \Tn\ and therefore / will always be greater or 
equal than both li and I2 so that (4.5a) will apply, and it gives q = I + 1 — Zq with 



[0 


for 


/ — m 


even 1 




for 


/ — \m\ 


odd j 



We also have have p = p' & Z"'"U{0} and 2a = 2f3 = \m\ and then the confluent hypergeomet- 
ric functions are just the generalized Laguerre polynomials: iFi{—p, \m\ + l,z) oc Lp"^^\z). 
Finally, because of the existence of the x ^ —x symmetry in the scalar case, we have that 
oBimiu = ^o^imuj in (3.23) and therefore the zero of the outer solution, if it exists, will be 
located exactly at x = 0. All these results for the scalar case coincide with [7], [8] and [9]. 
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V. NUMERICAL METHOD 

Two different methods have been used to obtain the numerical data. One method is the 
one used by Sasaki and Nakamura [17], consisting in approximating the differential equation 
(3.3) by a finite difference equation, and then finding the eigenvalue as the value of hEimuj 
that makes zero the determinant of the resulting (tri-diagonal) matricial equation. We have 
used this method to find the eigenvalues for several large values of c. However, we used the 
shooting method described in [18] to calculate the spin- weighted spheroidal function. 

The shooting method is applied in [18] to the spheroidal differential equation (i.e., /i = 0), 
and we adapted it to the spin-weighted spheroidal differential equation as follows. In general, 
for an initial, arbitrary value hEimu for the eigenvalue, which is different from the actual 
eigenvalue hEimuj, the numerically integrated solution is a combination of both the regular 
and the irregular solutions, i.e., 

where hUimw irregular solution at x = ±1, hUh^ is the numerically obtained solution 

and hUirm^ is the analytic, regular solution. A and B are unknown functions of hEimu- We 
need to modify the value of hEimu so that only the regular term Ahyimui is retained. In 
the scalar case, the boundary condition at X2 = +1 — dx may be imposed by requiring 
that h=oEimuj is a zero of the function g{h=oEimu;) = h=oy'iZT{^'i) - h=oy'imw{x2) , where the 
analytic value hy'imuji^'i) is known for the scalar case because h=oyimuii^) ^ h=oy[muj{~^) ■ The 
function g{h=oEimiu) should tend to zero as h=oEimuj approaches the correct eigenvalue and 
should tend to infinity when it is far from it because of the behaviour of the irregular solution. 
However, in general we have hy'imwi^) ^ ~-hy[mu}{~^) ■> relating solutions of equations with 
different helicity when h ^ 0, and therefore we do not know the analytic value hy'imuji^^) ^'^^ 
a particular value h ^ oi the helicity. We therefore decided to apply the shooting method 
by finding a zero of the function 

ginEim.) ^ ky'r-{x2) - .y;— (X2) (5.2) 

instead, where h.2/;mtj"^°^(^2) is not the actual analytic value, which we do not know, but an 
approximation to it: 

.y;rr^(^2) - ^f^.yL.(^2) (5.3) 
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Sasaki and Nakamura's method, which they only develop explicitly for the case h = —2 
and m = solves the angular differential equation (3.3) re-written with derivatives with 
respect to 6 rather than x: 

^ + [2{a -(3) + 2{a + /3) cos9 + cos9] ^+ 

(5.4) 

+ hEtmu -{a + p){a + p+l) + c^ cos^ 9 - 2hccos9 \hyimU9) = 

This equation is approximated by a finite-difference equation. Apart from at the boundaries, 
the derivatives are replaced with central differences. At the boundary points, the regularity 
condition (3.1) requires that dhyimuj/d9\^^^-^ = and the first order derivative (which has 
a factor 1/ sin^ in front) is approximated by a forward/backward difference atx = +l/ — 1 
respectively. The result is that equation (5.4) is approximated by 

{A9y 

+ [2{a -(3) + 2{a + (3) cos9, + cos 5,] 

sm fj ZlXo 

+ [hEimu. -{a + f3){a + f3 + l) + c^ cos^ 9, - 2hc cos 9i] hvL^ = 0, 

for i = 2, . . . , 2iV 

2(1 + 2a) + 

+ [h,^im^ -{a + f3){a + f3+l) + c^ - 2hc] nyLc. = 0, for z = 1 (^ = 0) 



(5.5) 



(Ae)2 

+ [hEimu. -{a + P){a + P+l) + c' + 2hc] nyLu, = 0, for ^ = 2iV + 1 (^ = tt) 

where 9i = 7r(? - 1)/(2A^) = A9{i - 1) and i = 1, 2, . . . , 2A^ + 1. Equation (5.5) can be 
represented as the product of a square, tridiagonal matrix A of dimension (2A^+1) x (2A^ + 1) 
and the vector of elements hVimuj equal to zero. In order to find the eigenvalue, Sasaki and 
Nakamura's method imposes that the determinant of matrix A is zero. 

We found that, already with = 100, for most modes the values of hEimuj obtained 
to quadruple precision actually provided values of the determinant so large that were even 
greater than the machine's largest number. We therefore decided to use this method only 
to find eigenvalues and use the shooting method when we wish to find both eigenvalues and 
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spherical functions. In fact, Sasaki and Nakamura's method without finding the spherical 
funcion is so much faster than the shooting method that the former is the preferable method 
to use if we wish to find eigenvalues far from any known eigenvalue (as we analytically do for 
c = for example). This is why we used Sasaki and Nakamura's method to find the eigen- 
values for large frequency and then used the resulting eigenvalue to find the corresponding 
spherical function with the shooting method. 

We wrote a program that implements Sasaki and Nakamura's method to find eigenvalues, 
particularly adapted to the case of large frequency. It calculates h^imoj rather than hEimu 
since hMmui ~ 0(c) for large c whereas hE imu) ~ ^l*-"^)' -^^ starts with the known value 
of h\,m,u)=Q (1-3) and finds the eigenvalue h^imuj for increasing frequency by looking for a 
zero of the determinant of the matrix A. This procedure is smooth no matter how large 
the frequency is if / < Zi or / < /2- However, if I > max(/i,/2), for some large value 
of the frequency, the eigenvalues for two consecutive values of / are so close (since they 
correspond to the same q and therefore their leading order term for large frequency is the 
same) that the initial bracketing of the eigenvalue includes both eigenvalues and therefore 
det A calculated with the values of h^imw at the two ends of the bracket has the same sign. 
From this value of the frequency on, instead of looking for a zero of the determinant the 
program just looks for the value h^imuj that is an extreme of the determinant. The reason is 
that this provides a point which is in between the two actual eigenvalues and it is therefore 
useful both as an approximation and as a bracket point for either of them. Instead of using 
minimization/maximization routines, which are very costly in terms of accuracy and time, 
in order to find an extreme of det A, the program looks for a zero of the derivative of det A, 
which can be calculated to be 

^^^ = trace [{det A) A-'] (5.6) 

and is very easy to evaluate. The program we have just described provided the graphs of 
h^imuj as a function of u for large frequency and there is therefore no need for it to distinguish 
with accuracy between the two consecutive eigenvalues. 

The extreme point of the determinant found by this program is used by another program 
which uses the shooting method and Runge-Kutta integration to bracket and determine the 
two close eigenvalues and their corresponding angular functions. It initially looks for a zero 
of the function g{hEim,uj) inside a bracket of the eigenvalue. If it finds a zero inside the 
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bracket, it then directly implements the shooting method as described in [18]. Instead, if 
it does not find a zero inside the bracket it then assumes that it is because the frequency 
is large enough so that there are two eigenvalues inside the bracket corresponding to two 
different, consecutive /'s. It then looks for a minimum of g{hEimS) (with a possible change 
of sign if there is a maximum instead) and uses that minimum to find a zero to its right 
or to its left depending on which one corresponds to the / we are interested in, according 
to (4.5). This second program also finds the zero of the function hSimu close to x = for 
large u if it has one as indicated by (4.6), uses a smaller stepsize in x close to x = ±1 to 
cater for the rapid oscillations of the angular function there for large u and makes use of 
equations (4.5) and (3.5) to help bracket the eigenvalue. This program provided the graphs 
of hSimA^{9) for large frequency. 

Both programs were written in FortranQO and contain parallel algorithms that use the 
Message-Passing Interface as the message-passing library. 

VI. NUMERICAL RESULTS 

All the numerical results and graphs in this section have been obtained setting a = 0.95 
and M = 1. 

There is an obvious numerical problem when p, p' G Z'^UjO}. In this mentioned in 

Section IV, the eigenvalues for two different values of / (but same h, m) become exponentially 
close as c increases ( [11]). This means that for this case we are not able to find the functions 
for very large values of the frequency. For example, in the case below for h = —1 and m = 1, 
when u = 25 the eigenvalues for / = 3 and / = 4 only differ in their 14th digit. 

BRW do give the analytical value for q for spin-0. For spin different from zero, however, 
they try to numerically match their large-frequency asymptotic expansion of the eigenvalue 
with the expansion for small frequency given by Press and Teukolsky ( [12] and [13]). As 
can be seen in Figures 2 and 3, this matching at intermediate values of the frequency might 
be good for certain cases, especially for small /, but not for other ones. All eigenvalues start 
off for frequency zero at the value given by (1.3), as expected, and when / > /i or I2 the 
pairs of curves that share the same value of q become exponentially closer and closer to each 
other as the frequency increases. When the frequency is as large as 100, the curves fully 
coincide in the expected pairs for large frequency (given by equation (3.5), and BRW for 
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lower order terms) where q comes in as a parameter. From this, the corresponding value 
of q for a certain set of values of h} can be inferred, and this coincides with the one 

given by equations (4.5a), (4.5b) and (4.5c). 

We calculated and plotted in Figure 4 the SWSH for /i = — 1, / = 3 & 4, m = 1, where 
the value of g, given by (4.5b), is the same for both of them: g = 4 (this is a case where 
p^p' G U {0}). Several features can be seen. Firstly, as the frequency increases from 

= 5 to 25, the functions become flattened out in the middle region of x and squeezed 
out towards the edges. Since the value of q is the same for both cases, the inner solution is 
the same for both of them, with the only exception of the relative sign between the inner 
solution for positive x and negative x (3.34). The function for / = 4 has three zeros and the 
one for / = 3 has two, in agreement with (4.1). The inner solution provides for the two zeros 
of / = 3 and the corresponding two of / = 4, and these become closer to the boundary point 
X = +1 as the frequency increases. The extra zero of / = 4 comes from the outer solution 
and becomes closer to x = with increasing frequency. 

In Figures 5-7 the lines labelled as 'inner' have been obtained with (3.9), the ones labelled 
'outer' with (3.23), the ones labelled 'uniform' with (3.38) and the ones labelled 'numerics' 
with the programs described in Section V. These figures show that the outer (normalized to 
agree with the numerical data at x = 0), inner (normalized to agree with the numerical data 
at x = ±0.96) and uniform (also normalized to agree with the numerical data at x = 0) 
solutions approximate the numerical data for u; = 25 in the boundary layers and in the 
neighbourhood of x = 0. The outer solution is valid until the boundary point x = — 1 but 
not until x = +1 since the function has two zeros close to it and the outer solution cannot 
cater for them, whereas the uniform solution is a valid approximation for all x. The inner 
solutions, on the other hand, prove to be a good approximation in the boundary layers but 
not close to x = 0. 

Figures 8 and 9 prove equation (3.34) to be correct for the case m > 1: for the specific 
values h = —1, / = 4, m = 1 and g = 4 the inner solution (3.9) has been normalized to 
match the numerical data at the points x = ±0.998 for different values of the frequency 
from 5 to 25, in order to be able to calculate _i-D4,i,a;, -iC'4,i,w and -i-D4,i,a;/-iC'4,i,a;- When 
plotting this numerical ratio together with the analytical result (3.34), the two lines are 
parallel and therefore agree to highest order, and the ratio between the numerical and the 
analytical data tends to 1. 
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Figures 10-15 correspond to modes with h = +2, m = 1, = 35 and / = 5 or / = 6. The 
modes for both values of / yield q = 6. However, the mode with / = 5 does not possess a 
zero at xq whereas the mode with / = 6 does. The behaviour for positive x is very similar 
for both values of / but for negative x the behaviours for the two modes differ by a sign. 

For h = —1, / = 2, m = 1 and uj = 100, the corresponding value of q is 2. This is a case 
where p G Z'^ U {0} and p' ^ Z"*" U {0}. The numerical solution together with the uniform 
expansion (3.32) is plotted over the whole range x G [—1, 1] in Figures 18-20. 

As we have seen, in this case the function has an exponential behaviour far from the 
boundary layers, so that a plot of the log of the function allows us to see the behaviour 
over the whole range of x. Both the uniform expansion and the outer solution have been 
normalized so that they coincide with the numerical value at x = 0, and the inner solution 
has been normalized once at x = 10~^ and once at x = —10"*^. The uniform expansion 
agrees with the numerical solution for all values of x. The outer solution agrees with the 
numerics everywhere except very close to x = ±1, where it veers off. The inner solutions are 
valid all the way from their respective boundary layers until, and past, x = 0, which is due 
to the exponential nature of the function in the region between the boundary layers. The 
inner solutions show a jump at x = due to the different orders in c of hCimuj and hDimu)- 

The above features can be seen in detail for x close to and ±1 in Figures 21-23 where 
they have been rescaled by 10^° for x close to and —1. 
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FIG. 2: _iA, function of uj for several / and q. The red crosses are the numerical data. 

The navy blue lines are using BRW's expansion for hMmuj and the light blue lines are Press and 
Teukolsky's. 
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FIG. 3: +2^l,i,ij as a function of uj for several / and q. The red crosses are the numerical data. 
The navy blue lines are using BRW's expansion for hMmuj and the light blue lines are Press and 
Teukolsky's. 
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FIG. 5: ~iSi^i^25 for Z = 3 & 4. Different solutions as labeled. The continuous lines correspond to 
I = 4 and the dotted ones to I = 3. 
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FIG. 6: -i5;^i^25 for / = 3 & 4. The curves above the x-axis correspond to / = 4 and below the 
axis to / = 3. Correspondence between colours and solutions is the same as in Figure 5. 
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FIG. 7: -i5';^i^25 foi^ / = 3 & 4. The continuous lines correspond to / = 4 and the dotted ones to 
/ = 3. Correspondence between colours and solutions is the same as in Figure 5. 
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FIG. 8: ^ c^'^''^ for w = 5 ^ 25. The analytic values have been obtained with (3.34). 
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FIG. 9: Ratio between numeric and analytic values of \c*^]^'^ ■ analytic values have been 

obtained with (3.34). 
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FIG. 10: +2'S'5,i,35- Green line corresponds to uniform solution (3.38) and red line to numerics. 
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FIG. 11: +2'S'6,i,35- Green line corresponds to uniform solution (3.38) and red line to numerics. 



35 



2.5- 
















2 














/ / jtWt, \\ 

//*t *\\ 
t / \t\l 


1.5- 














/ / / \*4 

i VI 

1 rl Vl 


1- 














/ ^/ li 


0.5- 






' .1 1 1. ir 













0.6 


0.65 


0.7 


0.75 


0.8 

X 


0.85 


0.9 0.95 1 



FIG. 12: +2'S'i,i,35 for / = 5 & 6. The continuous lines correspond to / = 6 and the dotted ones to 
^ = 5. Correspondence between colours and solutions is the same as in Figure 5. 
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FIG. 13: +2'S';^i^35 for / = 5 & 6. The continuous lines correspond to / = 6 and the dotted ones to 
/ = 5. Correspondence between colours and solutions is the same as in Figure 5. 
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FIG. 15: +2'S'6,i,35. Correspondence between colours and solutions is the same as in Figure 5. 
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FIG. 16: - 




/ = 5. The continuous lines correspond to the analytic expression (3.35) and the dotted ones to 
the numerical data. 
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FIG. 17: Ratio between numeric and analytic values of ^^^^''^'"^ for u; = 1 ^ 34. Blue lines (plotted 
both continuous and dotted to show agreement with red line for large tu) correspond to / = 6 and 
red line to / = 5. 
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FIG. 19: log(_i52,i,ioo)- The continuous, green line corresponds to the uniform solution (3.32) and 
the dotted, red one to the numerical data 
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FIG. 20: log(_i52.i,ioo)- The red line (numerical data) overlaps with the navy line (outer solution). 
The light blue line (inner solution valid at x ~ —1) and the magenta line (inner solution valid at 
X ~ -|-1) overlap with the red/navy lines for negative and positive x respectively. 
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FIG. 21: -i52, 1,100- 
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FIG. 22: 104"_i52,i,ioo 
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